Propagation of interacting force chains in the continuum limit 
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We study the effect of mergers in the force chain model describing the stress profile in static gran- 
ular materials. Combining numerical and analytical calculations we show that granular materials do 
not generally behave in an elastic-like manner, however they may under specific conditions, which 
are elaborated. Non-elastic behavior resulting from the non-linearity of the full force chain model 
is discussed. 
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A striking characteristic of stress transmission 
in granular matter is the network of highly singu- 
lar lines, termed force chains, along which stress 
propagates Q, |M 1^ Q - This force chain network 
reflects the specific packing of the system, which 
is unique for each experiment. Although this leads 
to significant fluctuations in stress profiles from ex- 
periment to experiment, the average stress profiles 
calculated over an ensemble of similar experiments 
seem well defined In particular, the ensemble 
average of the response to a small localized force, 
the response function j^, is a bell-shaped curve 
with scaling properties similar to that of an elastic 
response function 0;I3,0;I3- I* is natural, there- 
fore, to ask whether the average stress in granu- 
lar materials behaves according to elasticity the- 
ory. In this Communication we address this ques- 
tion within the context of the Force Chain Model 
(See for a discussion of a model of masses 
linked with linear and nonlinear springs.). We find 
significant deviations from elasticity, except for the 
case of an isotropic packing with nearly isotropic 
applied forces. 

The recently proposed Force Chain Model 
(FCM) transforms the singular behavior of 
stress in states of a granular material into a con- 
tinuum theory by averaging over an ensemble of 
states. This is done by writing a Master equation 
for the average density of force chains, allowing 
force chains to propagate, split and merge |^llol|. 
Previously, the response function of the FCM was 
calculated in three different fashions: simulating 
force chain propagation in small quenched disor- 
dered media calculating the constitutive rela- 
tion for granular materials on large scales using 
a splitting-only variant of the force chain model 
Q; and linearizing a specific discretized version 
of the model around a homogeneous solution 10]. 
The response functions calculated by the first two 
methods agreed qualitatively with the experimen- 
tally measured response 0, Q, exhibiting a bell 
shaped peak, while the third method gave a tran- 
.'sitinri frnm a. .'sinp'le tn a. double neak at an inter- 



mediate length scale [Tol |. 

A priori, there are reasons to expect granular 
materials to behave non-elastically: They cannot 
sustain tensile stresses, they rearrange when ex- 
ternal loads are changed, and they have no equi- 
librium stress-free state with respect to which to 
define a displacement field (See Q for a differ- 
ent view on these issues.). Despite these consid- 
erations, the central result of Reference ^] was 
that, in the absence of force chain mergers, gran- 
ular materials behave in a quasi-elastic manner on 
large length scales. In this Communication, we 
argue that the effect of force chain mergers is to 
change this: Generically, granular materials do not 
behave elastically. The apparent elastic-like be- 
havior found in experiments 0j 13 is restricted to 
specific packing geometries (i.e. isotropic) and to 
specific configurations of applied loads (i.e. near- 
isotropic) . 

Stress profiles in the FCM may be calculated 
in three ways, each adapted to a different length 
scale: Monte-Carlo simulation on small scales; nu- 
merical solution of the discretized model of Ref. 

on intermediate scales, and calculation of the 
constitutive relation of the full FCM on large 
scales. We will discuss the latter two methods; the 
simulation results will be presented in We 
will show that observed deviations from elasticity 
13 may be understood in the context of the FCM. 

In the framework of the FCM, a force chain is 
characterized by its intensity, / (the pressure ex- 
erted on each grain along the force chain), its direc- 
tion, n , (which is determined with respect to the 
applied force on the boundary), and its position, 
r. There are four events that involve the creation 
or annihilation a force chain {/, n} at r [iol |: it can 
split; another force chain can split, creating it as 
one of its offspring; it can merge with another and 
so be annihilated; or two other force chains can 
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for the force chain density, P = P{f, n, r) 0,^3- 



h-VP ^ ~-P 
A 



(1) 



+ 



^ j Pii^oSifh - ihhi + hn2))dhdh 
-QP J Pi^oKhn2 - ififii + fn))dhdh 
-| J PiP2^oS{fn - (Aril + hn2))dhdh. 



Here Pi — P{fi,ni,'r), Q is the force chain width 
[T^ . A is the sphtting mean free path of a force 
chain, and the functions ipo a-nd ipo are the weights 
of a sphtting/merging event (which depend, in 
principle, on the directions n, ni, and 71,2). dfj = 
dfjdhj^ and the delta functions ensure force bal- 
ance. Throughout this Communication we will as- 
sume that when two force chains meet, they merge, 
i.e. (po = l. 

Following [3 , we define the force chain intensity 
density, F{n,f) = J P{f,n,r)fdf and its angular 
moments: 



Ja{r) ~ J 'naF{n,r)d'h (2) 

o'Q/3(r) = aD j nanpF{fi,r)dh; (3) 

a is the grain size, and D the dimension of the sys- 
tem. (Jap is the local stress tensor 13] and J{r) can 
be thought of as the average force chain current. 

In order to gain insight on the mesoscopic 
scale of stress profiles in a granular material, 
we approximate Eq.QJ, following Ref. pH ]. 
We employ the discrete ordinate method pro- 
posed by Chandrasekhar 14] for solving the 
Radiative Transfer Equation. This approxi- 
mates the integrals in Eq. |^ by sums, by 
discretizing the directions of force chains as: 
P{f.fi,7^ = EtiP^Sif - r)<5(n-n,) where 
hi — {cos0i,sm9i),i = I, ..,6, and 0i^i — 9i = ^ 
(with 67 = 6*1). The P^'s are six different func- 
tions representing the weights of the force chains 
propagating in directions fii. Note that the choice 
6i+i — 9i = ^ implies that all forces have the same 
intensity, /*, in order to satisfy force balance. Sub- 
stituting this into Eq. |^ results in six coupled 
differential equations for the six force chain den- 
sities. Pi. Rescaling Pi jq^i ^^'^ r ^ ^ we 
arrive at the dimensionless equations |1C| : 



In Ref. 10] homogeneous solutions of the form 
{Pj} = {q, q^, q, q~^, q~^ , for any q were con- 
sidered; however there are others, for example: 
{Pj} = {q-\l,q,q,l,q-^}. In 10], Eq. was 
solved by linearization around the first homoge- 
neous solution, and, remarkably, a double peaked 
response was shown to emerge at intermediate 
depths. It is important to understand what phys- 
ical conditions these solutions correspond to. Let 
us begin by choosing q = 1 m. either of the ho- 
mogeneous solutions; this means that the force 
chain density, and therefore the stress, is uni- 
form throughout. This represents the case of an 
isotropic array stressed hydrostatically, which is 
the reference frame for response function experi- 
ments. In the case that g 7^ 1 we have some pre- 
determined relation between force chain densities 
in different directions at the same position which 
is unphysical for an isotropic material. 

Numerical solutions of the descretized model, 
Eq. Q), were obtained using a second-order ac- 
curate finite-difference approximation and solving 
the resulting nonlinear algebraic system of equa- 
tions iteratively ^3] • In Figure ^ the response 
to a normal force as calculated by this model is 
plotted 0. The peak is symmetric, as expected. 
In contrast, the response to a tilted force (Figure[21 
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FIG. 1: The response function as calculated with 
boundary conditions: Pi — 0. The width of the re- 
sponse scales linearly with depth as seen in the insert, 
where the response function at various depths is plot- 
ted; the curves normalized by peak height. See foot- 
note fl5| . 

shows the pressure profile on a plane normal to the 
applied force.) is asymmetric. The width at half 
maximum of both functions increases linearly, in 
accord with the experiments of Ref. '3] . The bell 
shaped peak of the response to the normal force 
is in agreement with the elastic-like behavior as- 
cribed to granular materials (see e.g. |3. ITiL 
However, the asymmetric response to the tilted 
force deviates from the elasticity prediction fisl ]. 



-VP, = -P,; + P,+i + 

+ Pi+lPi-l ~ Pi+2Pi — Pi-2Pi (4) 

Note that these equations are written for isotropic 
homogeneous assemblies, since the mean free path 
is assumed constant. 



An explanation for the deviation from elasticity 
of the response function to a tilted force can be 
found in Eq. ^ which connects the force chain 
density in one direction with that in any other di- 
rection. This means that if all force chains arriving 
at the surface are canceled (the grains rearrange to 




FIG. 2: The response of a granular assembly to a 
force tilted 60° with respect to the free surface. In- 
sert: Stress normalized by peak height, showing linear 
scaling with depth. See footnote [l^. 



have zero stress at a free surface, for instance) the 
total force chain density will be zero in the vicinity 
of that surface. For an elastic material, however, 
it is possible to have no strain in the direction per- 
pendicular to the surface and a finite strain paral- 
lel to the surface, since the strain components are 
independent. It is noteworthy that this deviation 
from elasticity is observed also in the splitting-only 
version of the force chain model (see [13 ) . 

One of the fundamental characteristics of the 
full force chain model is its non-linear nature (see 
Equations and Ql). In order to estimate the 
effect of this non-linearity we tested superposition 
by comparing the response to two different per- 
turbations, first applied simultaneously and then 
applied separately. Figure E^a) presents the re- 
sponse to two forces applied close to one another, 
and Figure|3Ib) two forces applied further apart. It 
is clear that while the effect of force chain interac- 
tion is significant in the former case, it is negligible 
in the latter. 
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FIG. 3: A comparison between the response to two 
forces applied simultaneously (dashed line), and ap- 
plied separately (solid line), a) the forces are applied 
close to each other, b) the forces are applied at a dis- 
tance. Lack of super pos ition obtains for (a) but not 
for (b). See footnote (l^. 



While these results are suggestive, the dis- 
cretized model gives incomplete understanding be- 
cause all forces in the system are equal. Thus, to 
better oiu: understanding of the effect of force chain 



mergers we calculated the constitutive relation on 
the macroscopic scale by calculating angular mo- 
ments of Eq. (P) , in the spirit of Ref . 8] . Multi- 
plying Eq. Q by /n and integrating over ]h we 
arrive at the force balance equation: 

V • cr = Fo 

where Fq is an external body force. In order to 
calculate the second moment of Eq. one has 
to compute the integral: 

/ = y PiP2fnani3ipoS{fn-{fini+f2n2))d£dfidf2 



It has been shown both in experiments |2Cl l2lL l22l | 

and in simulations 23] that the probability distri- 
bution of forces, P{f), has a maximum, and decays 
exponentially for larger forces. Thus, we approx- 
imate the above integral by considering small de- 
viations of fi and /2 from the intensity, fmax at 
which the P{f) is maximal. That is, we write: 



/l — fmax + Sfl 
/2 = fmax + Sf2 



(5) 



and neglect terms with high orders of Sfi in /. 
Moreover, we assume near-isotropy, by expanding 
F{n, f) in spherical harmonics and keeping only 
the terms 



aF(n, r ) 2± p -|- Dh ■ J 



D + 2 



-n ■ a ■ n 



(6) 



where a is the traceless part of the stress tensor 
and p = j^Tr{a} is the pressure. This gives a 
constitutive relation [Tl|: 

(Ta/3 = A B{(p)V ■ JSap + Jafi - C {fl) ■ J {r)5al3 

- D{{n)^J{r)p + {n)pJ{r)^) (7) 



where (n) = 



J nP{f,h,r)df_, Ja/s = 
and whose constants A, C, D 
are determined by the specifics of the granular 
packing. The function B{(j)) depends as well on 
the force chain density (/)(r) = / P{f,n,f)df [Tll |. 

For nearly homogeneous and isotropic systems, 
the terms which are products of J and (h) are 
smaller than terms linear in J. If they may be 
ignored, the constitutive relation reduces to: 



ac(3 = A B{(f>)V ■ JS^p + J, 



Q/3 



(8) 



This equation is formally equivalent to the con- 
stitutive relation of conventional elasticity p^ . 
Therefore, we can define two pseudo-elastic mod- 
uli: The pseudo-Poisson ratio, ly = j^g, and the 
pseudo- Young modulus, E ^ A{1 + i^). These de- 
pend not only on the geometry of the pile, but also 
on Dosition through the force chain density (t>(r). 



Therefore, the pseudo elastic behavior obtains only 
for nearly homogenous systems. 

Generally speaking, the constitutive relations 
calculated by the force chain model (Eq. are 
different from those of conventional elasticity; in 
particular, they are nonlinear. This nonlinearity is 
somewhat subtle, and holds for the ensemble av- 
eraged stresses. For a given packing, which does 
not change upon application of external forces, it is 
clear that superposition must hold, since the grain- 
scale equations of force balance are linear. Con- 
sider, however, the ensemble of stress states which 
are consistent with given set of external forces 
{^i}. We believe that this ensemble is statistically 
different from that ensemble which is compatible 
with a different set of external forces {^2}, even 
if {^1} and {^2} are very similar. Physically, this 
is reflected in the fragility of the material: re- 
arrangements occur when the external conditions 
are changed. Thus, the FCM predicts that if an 
ensemble averaged stress field is measured, then 
granular materials will exhibit nonlinearity in its 
response. 

In this Communication we have dealt mainly 
with the ensemble average of the stress profile in 
granular materials. However, the singularity of 
the force chains and the wide distribution of force 



chain intensity measured [ij suggest that it might 
be interesting to study force fluctuations, and the 
effect of frictioiL, in the framework of the FCM. As 
in References the existence of force chains, 

in the sense of a reasonably straight line of grains 
in contact, was assumed. It remains to be seen un- 
der what conditions this assumption is reasonable. 
We expect that for very hard grains (more pre- 
cisely, small stresses compared to the grain com- 
pressibility), force chains will exist. If this is the 
case, the effect of friction would be to change the 
details of the packing obtained, such as the scat- 
tering mean free path and persistence length, but 
not undermine the existence of force chains (In- 
deed, grains are frictional in all real experiments 
I^IMQ' Furthermore, friction might stabilize a 
force chain network, by allowing it to bear loads 
which would otherwise be "incompatible", this in 
turn can lead to a larger load regime for which the 
granular material responds linearly. 
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